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Abstract. We study analytically and numerically the generation of 
shock waves in a quasi one-dimensional Bose-Einstein condensate (BEG) 
made of dilute and ultracold alkali-metal atoms. For the BEG we use an 
equation of state based on a ID nonpolynomial Schrodinger equation 
(ID NPSE), which takes into account density modulations in the trans¬ 
verse direction and generalizes the familiar ID Gross-Pitaevskii equa¬ 
tion (ID GPE). Comparing ID NPSE with ID GPE we find quantita¬ 
tive differences in the dynamics of shock waves regarding the velocity 
of propagation, the time of formation of the shock, and the wavelength 
of after-shock dispersive ripples. 


1 Introduction 

In the past years the zero-temperature hydrodynamic equations of superfluids [T] have 
been successfully applied to investigate dilute and ultracold bosonic superfluids, i.e. 
Bose-Einstein condensates (BECs) made of confined alkali-metal atoms [2]. More re¬ 
cently it has been shown [3] that also fermionic superfluids in the BCS-BEC crossover 
can be modelled by hydrodynamic equations and their generalizations [1]. 

The formation of shock waves induced by an ad hoc external perturbation has 
been observed in BECs mm and in the unitary Fermi gas [S] . Motivated by these 
low temperature experiments on shock waves in atomic superfluids, some authors 
have theoretically analyzed the formation of shock waves in various configurations of 
bosonic |1()I11I12I13I14TT3] and fermionic superfluids |lbll7ll8| . In particular, shock 
waves of a strictly one-dimensional (ID) superfluid in the regime of quasi condensation 
have been studied in Refs. |10lll] by using the ID Gross-Pitaevskii equation (ID 
GPE). 

Here we extend their analysis by studying a cylindrically confined BEG, whose 
transverse width a is not frozen but depends on the longitudinal axial density p 
[19120) . and adopting the ID nonpolynomial Schrodinger equation (ID NPSE) [T3] . 
Indeed, ID NPSE is a more realistic model for a comparison with experiments, where 
quasi one-dimensional BECs are produced by imposing a very strong harmonic con¬ 
finement of characteristic length a_L along two directions. In this paper we find that 
the shock waves obtained by using ID GPE, based on the hypothesis of a strictly one¬ 
dimensional Bose-Einstein condensate, show quantitative differences with respect to 
the ones derived from ID NPSE. 
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2 Hydrodynamic equations for a quasi ID bosonic superfluid 


We consider a dilute bosonic superfluid at zero temperature that is freely propagating 
in the longitudinal axial direction z and is confined in the transverse directions {x, y) 
by a harmonic potential of frequency lvx and characteristic length a_L = {U/ 
with m the mass of a bosonic particle. Under the adiabatic hypothesis that the trans¬ 
verse profile of the superfluid has a Gaussian shape with a size depending on the longi¬ 
tudinal density , one finds the following dimensionless hydrodynamic equations 

for the axial dynamics of the superfluid 


dv 


dt 


+ |(H = 0 


dz 



= 0 , 


( 1 ) 

( 2 ) 


where p{z,t) is the superfluid axial density at time t, v{z,t) is the superfluid axial 
velocity at time t and /i(p) is the bulk chemical potential. Here length is in units of 
a_L, time in units of wj , energy in units of Tiuj±, and axial density in units of aj . 
Thus the axial dynamics of the confined superfluid depends crucially on the functional 
form of the bulk chemical potential, that is given by [50] 



(3) 


where 

= l + g p L'[-^] (4) 

pa‘‘ 

is the transverse width, g = 2aslax is the interaction strength with the repulsive 
(og > 0) s-wave scattering length of the inter-atomic potential and a_L is the harmonic 
length of the transverse external potential |2^. Here L[x] is the Lieb-Liniger function 
m, which is defined as the solution of a Fredhholm equation and is such that 


L[x] 


X — for a; <C 1 

2 

fora:>l 


(5) 


As discussed in Ref. m, under the condition g < 1, for gp <C g^ one gets the ID 
Tonks-Girardeau regime [55] where cr = 1 and 

M + 1 (6) 

with 1 the adimensional transverse energy. Instead, for gp ^ g^ one finds the (mean- 
field) BEG regime where tr = (1 -I- ypY^^ and 


= 




9P _ 1 

Jp 2 


^\/l + 5P + 


VI + 5P 


(7) 


In this BEG regime one can distinguish two sub-regimes: the ID quasi BEG regime 
for <C gp <C 1 where tr ~ 1 and p = gp + 1, and the 3D BEG regime for gp 1 
where cr = {gpY^'^ and p = V^PP E3- 

Notice that a different expression with respect to Eq. © has been heuristically 
proposed in Ref. [53]. However, the ID nonlinear Schrbdinger equation of Ref. [53], 
which is an alternative model with respect to NPSE to describe transverse effects. 
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gives a ID equation of state very close to the NPES one, Eq. and both are in 
very good agreement with the 3D equation of state of the full 3D Gross-Pitaevskii 
equation |19I23| . In Eq. ([7]) the term gpl\/\ + gp is the axial chemical potential while 
(l/2)(-^l + gp + 1/v^l + gp) is the transverse chemical potential, where (1 + gp)^!^ 
is the transverse width of the superfluid. 

Eqs. © and © with Eq. © can be derived [19124) from the 3D Gross-Pitaevskii 
equation [5] by using a variational wavefunction with a transverse Gaussian shape 
ESI- They are remarkably accurate in describing static and dynamical properties of 
Bose-Einstein condensates under transverse harmonic confinement. Usually in Eq. 
© there is also a quantum pressure (QP) term, given by Ij{2y/p)d'^^/pjdz^, and 
in this case the hydrodynamic equations are equivalent to a nonlinear Schrodinger 
equation, the so-called nonpolynomial Schrodinger equation (NPSE) [19]. NPSE has 
been used to successfully model cigar-shaped condensates by many experimental and 
theoretical groups (see for instance [2(1127) 1. It has been also used to study the role of 
a varying transverse width [^ in the Josephson effect of a Bose-Einstein condensate 
with double-well axial confinement [7^ . 

The QP term is negligible as long as the longitudinal width of the density profile 
p{z) is larger than the healing length ^(p). A simple estimation based on the compar¬ 
ison between QP term and interaction term shows that ^(p) = (1 -I- gp)^!^ j \J2gp. As 
previously discussed, Eq. © cannot describe the very-low-density regime of Tonks- 
Girardeau (op <C o^, with o < 1), where the system behaves as a ID gas of impene¬ 
trable Bosons [20121122) . 


3 Sound velocity and shock waves 

First we observe that the sound velocity Cs(p) can be obtained from the Eq. © by 
using the thermodynamics formula = p{dp,/dp) and is given by 


Ap) = 


/5P(1 


bp) 


(1 -b 0p)3/2 


( 8 ) 


see also [3D] . The sound velocity Cs (p) is the speed of propagation of a small (infinites¬ 
imal) perturbation of the initial condition of axially homogeneous system. By using 
Cs(p) the two hydrodynamic equations © and © can be written in a more compact 
form as 


dp 

dp 


(9) 

dt 



dv 

dv 

Cs(p)^ dp 

(10) 



p dx 


where dots denote time derivatives and primes space derivatives. 

As discussed by Landau and Lifshits [T], exact solutions of these equations can 
be found by imposing that the velocity v depends explicitly on the density p. In this 


way one has = = the hydrodynamic equations become 


dp dp dv dp 


( 11 ) 


dv dp dv dp Cs{p)^ dp 
dp dt ^ dp dx p dx 


= 0 


( 12 ) 
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We now impose that the two equations reduce to the same hyperbolic equation 


dp . .dp 

where 

/ ^ /X dv 

c[p) = v[p) + —p 

(13) 

(14) 

from Eq. (fTTI). but also 


C(rt=.(rt+ p , 

(15) 

from Eq. (IT^. It is quite easy to verify that, given a initial condition f{z) for the 
density profile, the time-dependent solution p(z, t) of the hyperbolic equation satisfies 
the following implicit equation: 

Piz,t) = f{z - c{p{z,t))t) , 

(16) 

which holds for regular solutions. In general, the initial wave packet f{z) splits into two 
pieces travelling in opposite directions depending on the sign of c{p). For simplicity 
in the following we consider the right-moving part only. 

To determine c{p) we observe that, from the equality of Eqs. (ITTl) and (fT^ we get 

dv Cs[pY fdv\ ^ 

dp^ P \dp) 

(17) 

from which 

dv _ Cs{p) 
dp P 

(18) 

and 

1 X r cs{p) 

v[p) = ~ dp , 

J Po P 

(19) 


where we impose that at infinity the initial density is constant, i.e. po = f{z = ±oo), 
and the initial velocity field is zero, i.e. v{po) = 0. Notice that in Eq. (IT^ we have 
chosen dv/dp > 0 on physical grounds and that the relation between the velocity and 
density in Eq. (18) is precisely Eq. (101.4) of Ref. [T]. This formula can be integrated 
by using Eq. dH). One finally has 

v{p) = <P{p) - 'P{po) , (20) 


where 

'1 13 3 

2’~2’4’2’ 4^ 

The function Fi [a, &i, & 2 , c, x, y] is the Appel hypergeometric function given by 


-Pip) = 2^Fi 


:9P, -9P 


( 21 ) 


CXD OO 

Ei[a,6i,&2,c,a;,y] = X! X! 

m—0 n—0 


{c)m+nm\n\ 


-x-y" 


( 22 ) 


where {a)m = r{a + m)/F{a) is the Pochhammer symbol with r(x) the Euler gamma 
function. The velocity c(p) follows directly from the velocity v{p) by using Eq. (11) 
and Eq. (14). It reads 

c{p) = v{p) + Cs{p) = P{p) - P{po) + Cs{p) . 


(23) 
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Fig. 1. Velocities v and c as a function of the scaled axial density gp, with grpo = 2. Solid 
lines: Eqs. and (I23l- Dashed lines: Eqs. (1^ and (Unj. 


Eqs. (PHI) and (1^ with ^(p) given by (1^ are the main results of the paper. Note 
that in the ID regime {gp <C 1), where Cs(p) = one finds 


v{p) = 2y/^-2y/^ , 

(24) 

and 


c(p) = ^y/^- 2 y/^ , 

(25) 

which are the results obtained by Damski [TT|. 



In Fig.[I]we plot the velocities v and c as a function of the scaled and shifted axial 
density g{p — p^), chosing gp^ = 2 to enhance the differences between ID and NPSE 
hydrodynamics. The solid lines are obtained by using the NPSE hydrodynamics, i.e. 
Eqs. (pUl) and (1^ . while the dashed lines are obtained by using the ID hydrody¬ 
namics, i.e. Eqs. (El and El- The two panels clearly show quantitative differerences 
between ID and NPSE hydrodynamics. Notice that v{po) = 0 while c(po) = Cs(po)- 


4 Formation of the shock wave front 

Up to now the initial shape of the wave has been arbitrary. We consider now an 
example by choosing the following initial density profile 

f{z) = pq + pop e-^"/( 2 <T)^ ^ ^26) 

where p describes the maximum impulse with respect to the density background po- 
Both amplitude Airj) and velocity 12 ( 77 ) of the impulse maximum are constant during 
time evolution. The amplitude is A{p) = po(l -b p) while the velocity reads 

V(77) = c(po(l + 77)) . (27) 

As expected, taking p = 0 the velocity of the impulse maximum reduces to the sound 
velocity c(po) = Cs(po). Moreover, bright perturbations (77 > 0) move faster than 
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Fig. 2. Scaled time of shock formation Ts /{2<j) as function of the scaled impulse maximum 
por]. a is the initial width of the impulse and we set gpo = 2. Solid line: Eq. dig. Dashed 
line: Eqs. ISSJ. 


dark ones [r] < 0) [TB]. Let us consider a bright perturbation. The speed of impulse 
maximum is bigger than the speed of its tails. As a result the impulse self-steepens in 
the direction of propagation so that the formation of a shock wave front takes place. 
The time Tg required for such a process can be estimated as follows: the shock wave 
front appears when the distance difference traveled by lower and upper impulse parts 
is equal to the impulse half-width 2cr, namely [V{r]) — V(0)]T;, = 2a. It gives 


Tg = 


2a 


c(po(l + ?7))-Cs(po) ’ 


(28) 


where the local velocity c(p) is given by Eq. while the sound velocity Cs(p) is 
given by Eq. ([8]). In the ID regime (gp <C 1) the formula of the time Tg reads 


Tg = 


2a 

3(V1 + »7 - 1)\/^ ■ 


(29) 


In the case of a dark perturbation (r\ < 0) the tails of the wave packet move faster than 
the impulse minimum and the time of shock formation is simply Tg = 2ajicg^po) — 

c(po(l +??))■ 

In Fig. [2]we plot the scaled time Ts/(2ct) of shock formation as a function of the 
scaled impulse maximum porj. Again we choose gpo = 2. The solid line is obtained 
by using the NPSE hydrodynamics, i.e. Eq. (EHH, while the dashed line is obtained 
by using the ID hydrodynamics, i.e. Eq. (1^ . Both curves give Tg/{2a) —>• -boo as 
gpov 0 . 

It is important to stress that our analytical results have been obtained without 
taking into account the quantum pressure term \l[2^')d‘^^jdz^ in Eq. (2). Initially 
the spatial scale of density variations is cr, which can be chosen greater than the bulk 
healing length ^(po) = (1 + !\/‘^9Po- As a result the QP term is negligible at 

the begining of time evolution. 


5 Numerical results and after-shock dynamics 

As the shock forms up density modulations occur on smaller and smaller length scales 
beeing finally of the order of the healing length. Damski m has numerically shown. 
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for the strictly ID Bose gas, that the effect of quantum pressure term is that of 
preserving the single-valuedness of the density profile by inducing density oscillations 
at the shock wave front. Thus, after the formation of the shock Eqs. m and @ are 
not reliable. To overcome this difficulty we include the dispersive quantum pressure 
term in the hydrodynamic equations, which become 


dv 


dz 



dp 

di 

1 


+ ip^) = 0 


dz 

52 


2yp dz"^ 


Vp = 0 . 


(30) 

(31) 


We stress that at zero temperature for a viscousless superfluid the simplest regular¬ 
ization process of the shock is a purely dispersive quantum gradient term, which is 
proportional to in dimensional units. Clearly, Eq. © is a first order equation while 
Eq. (pil) is not due to the quantum gradient term. 

Introducting the complex field '0(z,t) such that 

V^(z,t) = (32) 

and 

viz,t) = ^0{z,t) . (33) 

it is strightforward to show that Eqs. dsni) and (IHTl) are equivalent to the following 
one-dimensional nonlinear Schrodinger equation 



1 52 

' 2 ^ 


■Mi^r) 


ijj 


(34) 


Using the bulk chemical potential p,{p) of Eq. ([3]) with ([T]), this nonlinear Schrodinger 
equation is the time-dependent generalized Lieb-Liniger equation we introduced some 
years ago [20j to accurately describe an experiment on a Tonks-Girardeau gas of 
®^Rb atoms [32| . In the BEC regime, where Eq. (j3]) reduces to Eq. 0, the time- 
dependent generalized Lieb-liniger equation becomes the ID time-dependent non¬ 
polynomial Schrodinger equation (NPSE) 



1 92 ffIV'P 
■ 2 ^ v'l+sIV'P 


7 : \/l + 5lV'P + 





(35) 


which gives the familiar one-dimensional Gross-Pitaevskii equation 



1 92 

2 ^ 


■5l'0P + 1 V' 


(36) 


in the ID quasi BEC regime, where <C 1 and the transverse width a of 

Eq. (m becomes cr ~ 1 (i.e. a zz a± in dimensional units) [13] . 

We solve numerically Eq. (13611 by using a Crank-Nicolson finite-difference predictor- 
corrector algorithm [53] with the initial condition given by Eq. ([33]) and v{z, t = 0) = 
0. In fact, as also shown by Damski we have verified that the initial velocity field 
v{p{z,t = 0)) and v[z,t = 0) = 0 give practically the same time evolution. 

In Eig. [3] we plot the time evolution of shock waves obtained with po = 4, cr = 4, 
and rj = 0.2. The nonlinear strength is choosen as g = 0.5, such that gp = 2 and 
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Fig. 3. Time evolution of shock waves. Initial condition given by Eq. (I2nj with (7 = 4, 
rj = 0.2 and po = 4. The curves give the relative density profile p{z, t) vs 2 at subsequent 
times t. Solid lines obtained solving the time-dependent ID NPSE, Eq. (EHl), and dashed lines 
obtained solving the time-dependent ID GPE, Eq. (136 II . We choose the nonlinear strength 
g = 0.5. 


consequently the system is very far from the Tonks-Girardeau regime (characterized 
by gp <C g^). The figure displays the density profile p{z,t) = \ip{z,t)\‘^ at subsequent 
times. Note the splitting on the initial bright wave packet into two bright travelling 
waves moving in opposite directions. As previously discussed, there is a deformation 
of the two waves with the formation of a quasi horizontal shock-wave front. Eventu¬ 
ally, this front spreads into dispersive wave ripples. The figure shows that there is no 
qualitative difference between ID NPSE (solid lines) and ID GPE (dashed lines) in 
the physical manifestation of supersonic shock waves. Nevertheless, due to the dif¬ 
ferent equation of state, there are quantitative differences. Our numerical simulation 
confirms that the velocity of the maximum of the shock wave is larger for the ID 
GPE with respect to the ID NPSE. 


In Fig. 0] we plot a zoom of the density profile p{z, t) of the shock wave to better 
show the after-shock dispersive ripples. In the figure we compare the density profiles 
obtained with ID NPSE and 3D GPE by chosing t = 30 for ID GPE and t = 40 for 
ID NPSE: in this way the spatial position of the highest maxima of the two profiles is 
practically superimposed. Fig. 4 clearly shows that the wavelength of the dispersive 
ripples is larger in the case of the ID NPSE. We have verifed that this is a general 
feature by performing other numerical runs with different initial conditions. 
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Fig. 4. Zoom of the density profile p{z, t) vs z of the shock wave. Initial condition given by 
Eq. (I26II with ct = 4, ?7 = 0.2 and po = 4. Solid line obtained solving the time-dependent ID 
NPSE, Eq. (1361) at t = 40, and dashed line obtained solving the time-dependent ID GPE, 
Eq. (1361) at t = 30. We choose the nonlinear strength g — 0.5. 


6 Conclusions 


We have shown that the shock-wave dynamics (velocity and density ripples) in a 
quasi ID BEG, described by the ID nonpolynomial Schrodinger equation (NPSE), 
is quite different with respect to the shock-wave dynamics of a strictly-lD BEG, 
described by the ID Gross-Pitaevskii equation (GPE). The dynamics of the quasi 
ID BEG becomes becomes equivalent to the one of the strictly-lD BEG only under 
the condition grp <C 1, with g the interaction strength and p the one-dimensional 
axial density. In particular we have found that by using ID NPSE the maximum of 
the shock wave moves slower, the time of shock formation is longer, and dispersive 
ripples have larger wavelength with respect to ID GPE. For both ID NPSE and ID 
GPE we have obtained analytical and numerical solutions which could be compared 
with experimental data. Indeed, shock waves can be experimentally produced by 
engineering the initial density of the superfluid. Working with BEGs made of alkali- 
metal atoms and confined in the transverse direction by a blue-detuned axial laser 
beam, a choosen axial density profile can be obtained by using a blue-detuned (bright 
perturbation) or a red-detuned (dark perturbation) laser beam perpendicular to the 
longitudinal axial direction. In this way one can experimentally test the equation of 
state of the quasi one-dimensional Bose-Einstein condensate analyzing the dynamical 
properties of the generated shock waves. 
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